%%%%%%%%%%%%%%%%%
% Return factor model,
% Control for other variables proposed by referee
% Liuren Wu, liuren.wu@baruch.cuny.edu
%%%%%%%%%%%%%%%%%
nx=5+7;nxe=nx+2; nlag=12;
minN=200;winth=0.01; 

dds=Z.dds;
ids=Z.ids;
dd=Z.dd; T=length(dd);
IV=mean(Z.IVV(:,1:2,1),2); %ATMV
y=mean(Z.NRV(:,1:2,1,3),2)*100; %straddle return, daily updated hedged
HVV=Z.HVV; %historical vol at 3 horizons
x=[Z.Factors, Z.ControlVariables];

ccx=NaN(nx+2,nx+2,T);
Bfv=NaN(T,nx);Afv=NaN(T,1);
statv=NaN(11,nx+2,T);
sx=NaN(T,nx);bi=NaN(T,4);
R2v=NaN(T,1); VB=NaN(T,nx);

Bfev=NaN(T,nxe);Afev=NaN(T,1);R2ev=NaN(T,1);

parfor t=nlag+1:T
    [Bi,nobt,Af,Bf,R2s,stats,sxt,sxt2,bit,r2it,qxt,qyt,qvt,vbt] =FunDailyORSestimation(t,y,x,HVV,IV,ids,dds,dd,nlag,nx,winth,minN);
    
    Bfv(t,:)=Bf; R2v(t)=R2s;
    Afv(t)=Af;VB(t,:)=vbt;
end

fl=[Afv,Bfv,100*R2v];
Tf=sum(isfinite(fl));
flstats=[nanmean(fl)', nanstd(fl,0)',sqrt(Tf)'.*nanmean(fl)'./nanneweystd(fl)',sqrt(12).*nanmean(fl)'./nanstd(fl)', nanacf(fl,1)', skewness(fl,0)', kurtosis(fl,0)'-3] ;




disp('%Table 8: Control for other variables');

xnamel={'Intercept', 'Delta hedging cost','Stochastic volatility','Random jumps', ...
    'Stock reversal','Stock momentum','Size','Book-to-market','Option open interest','Option bid-ask spread','Stock Price', ...
    'Historical premium','Volatility premium','Adjusted $R^2$'}';

fprintf(1, '  \\\\ \n');
for j=[1:4,nx,nx+1];fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,7),' \\\\ \n'], flstats(j,:));end
fprintf(1, '  \\\\ \n');
for j=[5:8];fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,7),' \\\\ \n'], flstats(j,:));end
fprintf(1, '  \\\\ \n');
for j=[9:11];fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,7),' \\\\ \n'], flstats(j,:));end
fprintf(1, '  \\\\ \n');
for j=[nx+2];fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,7),' \\\\ \n'], flstats(j,:));end

return
